A FINITE DIFFERENCE METHOD FOR THE 
WIDE-ANGLE 'PARABOLIC EQUATION IN A WAVEGUIDE 
WITH DOWNSLOPING BOTTOM 

D.C. ANTONOPOULOUtH, V.A. DOUGALIS*^ AND G.E. ZOURARIS^H 

Abstract. We consider the third-order wide-angle 'parabohc' equation of underwater acoustics in a cyhn- 
drically symmetric fluid medium over a bottom of range-dependent bathymetry. It is known that the 
initial-boundary-value problem for this equation may not be well posed in the case of (smooth) bottom 
profiles of arbitrary shape if it is just posed e.g. with a homogeneous Dirichlet bottom boundary condition. 
In this paper we concentrate on downsloping bottom profiles and propose an additional boundary condition 
that yields a well posed problem, in fact making it L'^-conservative in the case of appropriate real parame- 
ters. We solve the problem numerically by a Crank-Nicolson-type finite difference scheme, which is proved to 
be unconditionally stable and second-order accurate, and simulates accurately realistic underwater acoustic 
problems. 



1. Introduction 

We consider the third-order wide-angle 'parabolic' equation of underwater acoustics in a cylindrically 
symmetric fluid medium, [71 El [171 E] 

(WA) {l + qf3)vr + a'^q v,,r = i A (a^ +(3v), 

posed for {r, z) G 2? := {{r, z) G R^, < z < s{r), < r < R} for a given R > and a given bottom 
profile s — s(r). Here v — v{r, z) is a complex- valued function of the range r and the depth z, representing 




Figure 1. Domain of vahdity of (WA). 



the acoustic field generated in the fluid medium ('water') I? by a point time-harmonic source of frequency f 
placed on the z-axis, cf. Figure [H In (WA) we have put a = where /co — is a reference wave number 
and Co a reference sound speed, and A = where p, q are complex constants. The function f3 — (3{r, z) is 
complex- valued and represents — 1, where n is the index of refraction of the medium. 
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In practice /3 is real- or complex-valued with a small nonnegative imaginary part modeling attenuation 
in the water column. The coefBcients p, q are such that the rational function is an approximation 

to \^l + x near a; = 0. The p.d.e. (WA) is obtained formally as a corresponding paraxial approximation 
of the outgoing pseudodifferential factor of the Helmholtz equation written in cylindrical coordinates in 
the presence of azimuthal symmetry, [7], [T7]. The choice p — i, q — yields the standard, narrow-angle 
Parabolic Equation, [18], while the (l,l)-Pade approximant of y/1 + x, given by p = |, q = |, gives the 
Claerbout equation, |9]. In general, we shall take p and q complex; the choice p = 9+ |, Iin(q) < 0, [10^, has 
certain theoretical and numerical advantages as will be seen in the sequel. Let us also note that although 
in this paper we have in mind the application of the (WA) in underwater acoustics examples, our analytical 
and numerical methods can also be applied to wide-angle seismic wave, OE], and aeroacoustic f8| wave 
propagation. 

The p.d.e. (WA) is posed as an evolution equation with respect to the time-like variable r and the space- 
like variable z in the domain 2?, that has a range-dependent bottom described by z = s{r), where s is a 
smooth positive function on [0,i?]. We shall supplement (WA) by an initial condition modelling the sound 
source at r = 0, i.e. require that 

(1.1) v{0,z) = vo{z), 0<z<s(0), 

where vq is a given complex- valued function defined in [0, s(0)], and by the homogeneous Dirichlet boundary 
conditions 

(1.2) u(r,0)=0, 0<r<i?, 



(1.3) v{r,s{r)) =0, < r < i?, 

corresponding to a pressure-release surface and an acoustically soft bottom, respectively. 

There is considerable theoretical and numerical evidence to the effect that the initial-boundary- value 
problem (ibvp) consisting of (WA), and (|l.ip - (|1.3l) is well posed if the bottom is horizontal or upsloping, 
i.e. when s(r) < in [0, i?], and that it may be ill posed if the bottom is downsloping, i.e. if s(r) > for 
r G [0, i?], [14], [11], [12], [5]. In Refs. [Tl] and [12] an additional bottom boundary condition was proposed, 
that together with (WA) and (|l.ip - (|1.3|) yields, under certain hypotheses, a well posed problem for any 
smooth profile s. 

In [S] the authors of the paper at hand, in collaboration with F. Sturm, presented other types of additional 
bottom boundary conditions that render the resulting ibvp well posed and in addition, for real (3 and q, i^- 
conservative, in the sense that 

(1.4) / \vir,z)\^dz^ / \voiz)\^dz 

Jo Jo 

holds ioi < r < R. Specifically, it was observed that the ibvp consisting of (WA) and (ll.ll) - (|1.3p is 
i^-conservative, for real /? and q, if and only if the following boundary value condition holds 



0, 0<r<R, 



(1.5) Im J^{vz;r,s{r)) J^{v]r,s{r)) 
where for (r, z) e I? 

(1.6) J'iv, r, z) := q Vr{r, z) — i A v{r, z), 

provided that v satisfies (WA) and (|l.l[) - (|1.3p . Our main motivation for studying boundary conditions for 
which the 'energy' integral Jq*'"'^' |i'(?', z)P dz is conserved (for real (3 and q) in the case of the noncylindricral 
domain T) is the fact that this happens for solutions of (WA), (|l.l[) - (|1.3p in the case of a horizontal bottom 
and also for the standard PE {p — ^, q — 0) posed on T) as an ibvp with (ll.ll) - (|1.3p . for a general profile s 
and real (3. 

Here, we consider the ibvp (WA), (|l.ip - (|1.3l) for a downsloping bottom, with the additional boundary 
condition 

(1.7) vr(r,s{r)) =Q, < r < R, 
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which is equivalent to the condition Vz{r,s{r)) = 0, r G [0, i?], (in the case of a difFerentiable bottom with 
s(r) 0), as it seen by differentiating both sides of (11.31) with respect to r. Obviously, (ll.7p . in the presence 
of (11.31) . is also equivalent to J-{v;r,s{r)) = 0, i.e. satisfies (jl.Sp and yields an i^-conservative problem for 
/3, q real. In section 2 of the present paper, we prove that the resulting ibvp consisting of (WA), (ll.ip - (|1.3p . 
(jl.7p is stable in L^, and H^] also we show an _ff ^-stability result for the solution of the ibvp (WA), 
(jl.ip - (|1.3|) with an upsloping bottom. 

In order to develop a numerical method for ibvp (WA), (|l.l|) - (jl.3|) . (|1.7|) when the bottom is downsloping 
and g 7^ 0, we transform it, using the range-dependent change of the depth variable y = that renders the 
bottom horizontal, into an equivalent problem on the strip 0<?/<l, 0<r<i?. Letting := [0, i?] x [0, 1], 
I := [0,1] and u{r,y) := v{r,ys{r)) for (r, y) G fl, the ibvp consisting of (WA), (ITTT|) - (fr3| . (fTTT)) takes the 
following form 

A(r) (^Ur -i^u-^yuy'j = i u V (r, y) G n, 

(18) w(r,0)=0 Vre[0,i?], 

u{r,l) =Uy{r,l) ^0 Vre[0,i?], 
u{0,y) = UQ{y) ■.= vo{ys{0)) Vyel, 

where ^{r,y) :— 6(r, ys(r)) for {r,y) e and, for r G [0,i?], A(r) : H'^{I) n Hq{I) -> L^(/) is an indefinite 
one-dimensional elliptic differential operator in the y— variable defined by 

(1.9) A(r)«:=-t;"-ii±^^^^i^i-, V v e H^il) D H^{I), 

keeping in mind that the term — i u — y Uj, vanish at the endpoint of /. In section 3 we provide 
some conditions on the data of the problem that ensure invertibility of A(r) for r G [0, R] along with some 
regularity properties. We note that the p.d.e. in (|1.8|) . which follows from (WA) after the aforementioned 
change of variable, is not a usual Sobolev-type equation (cf. [TB]) like (WA) over a horizontal bottom (cf. 
[3], [1]). Due to the presence of the term ^ y Uy, the differential equation is of third order with respect to 
the space variable y and this offers an explanation of why an additional boundary condition may be needed. 

In section 4, for the approximation to the problem (II. 8p we propose and analyze a numerical method that 
combines Crank-Nicolson time-stepping with a standard second-order finite difference method in space. We 
would like to stress that the convergence analysis of the proposed finite difference scheme is not a repetition 
of the corresponding analysis for the flat bottom case (cf. [1], [3]). This is due to the fact that the differential 
operator is third order with respect to y, which leads to a truncation error of 0(1) at the nodes adjacent to 
the endpoints of /. Building up a careful consistency argument is important in proving a second-order error 
estimate. 

Finally, in section 5, we verify the accuracy and stability of the finite difference scheme by means of 
numerical experiments and apply it to solve an underwater acoustics problem in a downsloping benchmark 
domain comparing the results with those obtained from the model proposed in jl2l . 

For rigorous error estimates on other finite difference and finite element methods for the PE and WA 
equation on domains with horizontal or variable bottom we refer the reader e.g. to the papers [1], [3], [4], 
[6], [2], and their references. 

2. A PRIORI ESTIMATES 

Our aim in this section is to establish some a priori estimates for the solution v of the wide-angle equation 
(WA) on V under some or all of the auxiliary conditions considered in the previous section, and under several 
hypotheses on the coefficients and the bottom. In what follows we shall assume that the functions v, (3, s, and 
vq are sufficiently smooth so that the various estimates are valid. Also, in our analysis we shall employ L^, 
and range- dependent norms on [0,s(r)] which are given, respectively, by := (/q**''^^ \v{r, z)\'^dz)i , 
\\v\\i:=i\\vV + \\v4')'^ and H^l^ := (||^||2 + ||^,,||2)|. 

We begin by establishing some basic identities. 
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Lemma 2.1. Let q ^ and T he defined by p.6p . If v satisfies (WA), (jl.ip and (|1.2p . then, for <r < R, 
the following identities hold: 

f ||^;(r,.)|P=^(r)|«(r,.(r))p + fIm(g)||^;.(r,.)|P 

(2.1) 



(2.2) 



(2.3) 



s(r) 



Im [/3(r, z)] iTiv; r, z)\^ dz - ^ Im F{vz] r, s(r)) Tiv; r, s(r)) 







+ 2 Re 



i+gffl^'.z) ^^(i-^ z)ti(r, z) dz 



2AIm(i)||«,(r,.)|p 



i J'(i;^;r,s(r))i;(r,s(r)) 



-||«,.(r,.)f =i(r)k..(r,s(r))p + ^Re 



1 
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+ 4A 

(r) 



Im 



Vrir, ■)v^zir,z) dz 



2 Aim 



(I) Ik.. 



/3(r, z) T{v; r, z) Vzz{r, z) dz 



(1 + g/3(r,z)) |w.(r,z)|2 rfz ^qa^ |li;.^(r, OH 



(2.4) 



+ ia^ A 

- qa^ Vrzir,s{r))vr{r,s{r)) 



^v{r,z) + Vzz{r,z) Vr{r,z)dz 



Proof. We multiply equation (WA) using ()1.6p by J^{v;r,-), integrate by parts with respect to z and take 
imaginary parts to get (|2.ip . We then multiply (WA) using (|1.6p by w(r, •), fzz('', •)> respectively, integrate by 
parts, and take real parts, obtaining (|2.2p . (|2.3p . respectively. The last equality p.4p follows, by multiplying 
(WA) using (|1.6p by 'i;r('', •) and integrating. □ 

From p.ip it follows that if v satisfies (WA), (jl.ip . (|1.2p . and (jl.Sp . and /3 and q are real, then v satisfies 
the conservation property (|1.4p if and only if (jl.Sp holds. 

Before embarking on our study of the downsloping bottom problem, we prove a iJ^-stability result in 
the upsloping bottom case assuming only the homogeneous Dirichlet boundary condition (11.31) . thus com- 



plementing the i/^-estimate of [12]. The proof requires that i max s(r) be sufficiently small, i.e. a 'small 

frequency-shallow water' assumption. In what follows, c or C will denote generic positive constants, not 
necessarily having the same values in any two places. 

Proposition 2.1. Suppose that s{r) < for r G [0,i?], q ^ 0, and — max s{r) is sufficiently small. Then, 
there exists a constant c such that any solution v of the ibvp (WA), (|l.ip - (|1.3p satisfies 

(2.5) lk(r,-)l|2 <c||z;o||2, < r < R. 
Proof. We first prove that for re [0,R] we have 

(2.6) \\vr{r,-)\\+mv;r,-)\\< c\\v{r,-)h. 
To see this, note that gives for < r < i? 

(2.7) \\vrz{r,-W< ^\\vrirr)r + c\\vArr)\\\\vir,-)h + \vrzir,.s{r))\\vrir,s{r))\. 
Since Vr(r,0) — 0, we have 

(2.8) \\vrir,-)\\< sir)\\vrzir,-)\\, < r < R. 

By differentiating with respect to r the Dirichlet boundary condition v{r, s{r)) = 0, we obtain by the trace 
inequality 

\vr{r,s{r))\ = |s| \vz{r,s{r))\ < c ||w(r, •)||2. 
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Therefore, (P?7)) gives for < r < i? 

(2.9) Wvr^r, .)f < ^ ||«.(r, ■)r + c ||«„(r, .)|| ||z;(r, .)||2 + c |«..(r, s(0)l lk(r-, Olb- 

The equation (WA) solved for 

Vrzz B.nd ()2.8|1 yields now 

\\vrUr,-)\\ <c(||«.(r,.)|| + |k(r,-)ll2) 
<c(||«„(r,.)|| + ||i;(r,.)||2). 

Thus by Sobolev's inequality we obtain 

\vrz{r, s{r))\ < c ( \\vrz{r, ■)\\ + II 

<ci\\vr,ir,-)\\+\\vir,-)\\2). 
Using this in (12.91) and using (|2.8p we obtain, for any e > and < r < R 

WvrAr, Of < ^ s'ir) \\vr.{r, Of + e ||«..(r, .)|p + ||z;(r, .)||^. 
Hence, if — max s(r) is sufficiently small, we see, using again (|2.8p . that ||wr('^, Oil — ^ ||i'(r, OII2 for < J' < 

" 0<r<i?, 



i?. Hence, p.6p follows, since 

||^(^^;r,OII <c(||^;.(r,.)|| + ||^(r,0||) 

<c||'y(7-, OII2. 

Adding now (|2.2p and (|2.3p and using the fact that s(r) < 0, estimates from the proof of (|2.6p . and the 
Poincare-Friedrichs inequality we obtain for < r < R that 

|-(||«,(r,Of + ||«..(r,Of ) <c(||t;.(r,Of + ||«.(r,Of + ||«..(r,Of + ll-F(«;r,Of + ll«(^,Of ) 

<c(||«,(r,Of + lk..(r,Of ). 
Hence, from Gronwall's inequality we obtain 

\Mr, Of + |k..(r. Of < c(||i;.(0, Of + Of) , 

from which our conclusion follows. □ 

The 7?^-stability estimate (j2.5p implies of course the uniqueness of solution of the ibvp (WA), (|l.ip - (|1.3l) 
under the hypotheses of Proposition 12. 11 and forms the basis for a well-posedness theory for this problem. 

We turn now to the downsloping bottom case, with which we shall be concerned for the sequel of this 
paper. We first establish the following basic a priori estimates. 

Theorem 2.2. Supose that s(r) > for r E [0, i?], q ^ 0, and that either lm{q~^ + [3{r,z)) is of one 
sign in D or that ^ max s{r) is sufficiently small. Then the ibvp (WA), (jl.ip - (ll.3p . (|1.7p is L^-, H^-, and 

-stable. 

Proof Let v be the solution of the ibvp (WA), prT|) - ([0)) . pTTl) . Then it follows from ([TT^ that f{r,z) := 
F{v] r, z) satisfies 

/(r,0) = /(r,s(r))=0, < r < R. 

For each r e [0, R] consider the operator C{r) : n defined by 

C{r)u := + /3(r, 0) u + dfu, Vu G H^nH^. 

Then, under our hypotheses, £{r) is invcrtible in D Hq, in the sense that C{r)u — implies it = 0. To 
see this, note that from Cu = for u G H iJg, it follows that j^'"^^ Cuudz — 0, and by integration by 
parts, that 



(2.11) / {q-^ + P)\u\^ - a^\u.\^ 



■'(r) r 



dz = 0. 



Taking real parts in (|2.1ip and using the fact that < s(r)||M2||, it foUows that 
\\u,{r,-)f (l-X max s^{r) max |Re((j-i + /3)| ) < 0. 

\ 0<r<i? z^T} J 

Hence, if ^ max s{r) is sufficiently small we have that u = 0. Alternatively, taking imaginary parts in p. lip 

we obtain that j^'^^^ Im(g^^ + (3)\u\'^dz = 0, from which if Im(q^^ + /3) is nonzero and of one sign in I?, then 
u = Q. 

We conclude, from standard elliptic p.d.e. theory and the Fredholm alternative, [13], that given w ^ L^, 
then C^^w G n and ||£^^u'||2 < C||w||, for some positive C = C(r), independent of w. Since the 
coefficients of £ are smooth, C may be taken as a continuous function on [0, R\. Applying this result to the 
bvp ([TTU)) we see that for < r < i?, 

(2.12) ||/(r,-)l|2< C||z;(r,.)||. 
By the definition of /, it follows that on [0, i?] 

(2.13) ||«.(r,.)||, < C||i;(r,.)||,, £ = 0,1,2. 

Hence, by ()2.ip . L^-stability of the ibvp under consideration follows. To get the i?""^ -estimate, we use (|2.2p . 
()2.13p and the fact that Vz{z, s{r)) = for < r < i?, as remarked in the Introduction. It follows that for 
< r < i? 

i^\\v.{r,-W < C {\\v{r,-W + \\v.{r,-W) , 

from which _ff ^-stability follows by the Poincare-Friedrichs inequality. Finally, to get iJ^-stability, note that 
(1231) . ^^rm . (l2A3ll yield for < r < i? 



(2.14) f < s{r)\v..{r, s{r))\^ + C {Mr, Of + hz.{r, -W) . 

From Vz{r,s{r)) = for < r < R, it follows that Vzrir, s{r)) + s(r)vzz{r, s{r)) — 0. Since, Vzri'i'jsir)) — 
i /(r, s(r)), it follows by the trace inequality and (|2.12l) that for < r < i? 

^(r)k..(r,.(r))p=pi_K(r,.(r))p 

<C\\f{r.-)\\l 

<C\\v{r,-W. 

From (|2.14p the iJ^-stability estimate follows now. □ 
Remark 2.1. // Im(/3) > and lTii[q) < or if Im(/3) > and q is a real, nonzero constant, then 



Im (^^ + /3j > follows. The ibvp (WA), p.ip - (|1.3p . (|1.7p is of course L'^ -conservative if q and /? are real. 

3. Invertibility conditions for A 

Assuming that the operator A(r) defined by (|1.9p is invertible for all r £ [0,R] and that uq G Hq{I), we 
may write the problem ()1.8p equivalently as: 

Ur- ^yuy=i^^2p- T{r)u + i^u V(r,j/)efi, 
(3.1) M(r,l) = Vre[0,i?], 

u{0,y) = uo{y) VyEl, 
where T{r) := A"^(r) for r e [0,i?]. To simplify the notation we set 

"V i — ^y^^y) — ^ — ' ti'^j .- -^T^ 

for r e [0,i?] and y e I. 

Remark 3.1. Let g{r,y) := Ur{r,y) — 5{r)yuy{r,y). Assuming that the solution of p.ip is smooth on fi, 
we obtain the compatibility conditions g{r, 1) = g{r, 0) = 0, which yield that Ur(r, 0) = and Uy{r, 1) = for 
r € [0,R\. Finally, since u(0, 0) = uo(0) = 0, we get the surface pressure release condition u(r, 0) = for 
r G [0,R]. 
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Remark 3.2. Whereas the ibvp (WA), (jl.ll) - (jl.3l) . (jl.7p in the r, z variables is -conservative, in the sense 
that its solution satisfies (11.41) for q, j3 real, the solution of the transformed ibvp (jl.Sp in the r, y variables 
conserves the quantity s{r)\\u(r, ^ for q, 7 real. 

In the rest of this section we present some conditions on the data that ensure invertibihty of the operator 
A. The conditions are similar in nature to those in the hypothesis of Theorem l2.21 but it is useful to present 
them here because they are expressed in terms of the y variable and motivate analogous sufficient conditions 
for the invertibihty of the discrete operators in the next section. We use the notation j > 0, for the 

norm in the Sobolev space H^{I) and put = |jwj,||o,i. 



Lemma 3.1. Assume that 

(3.2) C^u inf 



[Cp 



^ (Re(g) + |g|2Re(7)) 



> 0, 



where Cpp > is the constant of the Poincare-Friedrichs inequality on I, i.e., \\v\\q 1 < Cpp |w|i,f for 
V € Hq{I), or that there exists S-i, G {1,-1} such that 



inf 



(3.3) C^s : 

Then, there exists a constant C > such that 
(3.4) 



(Im(g)-|q|2lm(7)) 



> 0. 



max ||T(r)V'||i,. < Cll^llo.z, y^beL^{I). 

re[0,R] 



Proof Let r € [0,R], v G H^il) n H^il) and e L^{I)- First, observe that 



(3.5) 



Re{A{r)v,v)oj ^ \\v' 



dy. 



When p.2|) holds, use of the Poincare-Friedrichs inequality and (|3.5|) gives 



(3.6) 



Re{A{r)v,v)o.i > 

>Cps\\v 
When (13. 3p holds, we have 



{CpfY 
2 



Rc(q+|g|S)^"('-) 



\v\Uy 



Im(A(r)u, u)o,j I >Si,lin{A{r)v,v)Qj 



(3.7) 



Im(9+|g|"-)s''(r) 



dx 



Since, by definition, T(r)?/' = A-'(r)-0 € H'^{I) n (/), we set v = T{r)ip in or dSJ]) and apply the 
Cauchy-Schwarz inequality to conclude for r £ [0, R] 

(3.8) \\T{rMoj < C\m,,,. 

Then, we combine (|3.5p for v ~ T{r)tp, and p.8p to get 

\T{r)^\l^ =Rei^,T{r)^),^, ^ 

^^■^^ <\mo^A\T{rm,.s + C\\T{r)^\\l 

<c\mi. 

Now dSH) follows easily from dSJ]) and (IXS)) . 



dy 



□ 

Assuming that (|3.2p or p.3p holds and using a induction argument, it is easy to establish that for m e Nq, 
there exists Cm > such that 



(3.10) 



max |r(r)/U+2,oo,. < Cm\f\m,^,i, V/e C"(/;C), 

r£[0,R] 



where on C™, \f\m oo i = inax|/'-™'(y)|. In addition, for to G Nq and £ G N, there exists d m > such that 

' ' yei 

(3.11) max |a^(r(r)/)|„,+2,oo,. < C,,,„ |/Uax{™-2.o},oo,., V/ G C'"-''{'"-2,o} 

rG[0,iJ] 

Remark 3.3. Differentiating both sides of the p.d.e. in p.ip once with respect to y, taking y ^ 1 and using 
the boundary conditions we get 

uyy{r,l) = -i^^^prl-^ {Tir)u)y{r,l), VrG [0,i?]. 

This reminds us in some sense of the boundary condition proposed in [12 for downsloping bottoms, which 
has the form 

Uyy{r, 1) = i I s(r) s(r) Uy{r, 1), VrG [0, R]. 

Remark 3.4. Condition p.2p follows from the hypothesis that ^ is sufficiently small for < r < R. If 
Im(7) > in Q and lTa{q) < or Im(7) > m 51 and q is real and nonzero, condition p.3p is valid for 
5* = -1. 

4. A Finite difference method 

In this section we construct and analyze a finite difference method for approximating the solution of the 
ibvp (IXTj) . 

4.1. Notation and preliminaries. Let J G N with J > 3. We introduce a partition of [0, 1] with width 
h := J and nodes yj := j h for j — 0, . . . , J. Taking into account the homogeneous Dirichlet boundary 
conditions at the endpoints of / we define the space of the finite difference approximations by 

On Xh we define a discrete L'^{I) norm || • ||o,h given by ||f||o,/i := {^^'jZri [''^jp) for ^ which 
is produced by the inner product {■r)o,h defined by (w,u))o,/i h'Y^''-Z^-^VjWj for v, w G Xh- Also, on 



Xh we define a discrete H^{I) seminorm | • \i^h by \v\i^h '■— ( hj^'^^ 



Vj + l-Vj 

h 



2\ 2 



for V G Xh, a discrete 



norm |1 • \\i^h by := ll«llo,/i + l^li.ft ' f^^^' ''^ ^ ^^d a discrete L°°{I) norm | • |oo,/i by 

|w|oo,/i := maxi<j<j^i for v G X/j. 

We define further the second-order difference operator A^, : Xh — > Xh by 

{Ahv), ^ j = l,...,J-l, V^GX,, 

and for r G [0, R], the discrete elliptic operator A/i(r) ; Xh Xh by 

{Ah{r)v)r.^~{^hv),-s^r)'-±^^^v„ j = l,...,J-l, Vz; G X^. 

Also, we define the first-order difference operator dh ■ Xh — > Xh by 

(dhv), := ^-'^'^h'-' , J = 1,...,J-1, VveXh, 

and introduce the auxiliary operators Ih : AT^ — > A;i given by {Ihv)j :~ for j = 1, . . . , J — 1 and 

V G A/i, and (g) : A^ — > A/i defined by {v (E) w)j — Vj Wj for j = 1, . . . , J — 1 and v, w E Xh- Also, we let 
a; G Aft be such that cUj := yj for j = 0, . . . , J — 1. For / : [0, 1] — > C wc define Phf G Xh by (Phf)] fiVj) 
for j = 1,..., J - 1. 
It is easy to show that 

(4.1) \\v\\o,h < "^Mi^h yveXh, 

(4.2) \v\oo.h< \v\i,h yveXh, 

(4.3) \v\^,h < h-^ \\v\\o,h V« G Xh, 
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and 

(4.4) {AhV,v)o,h = -\v\lh VweX,,. 

Let N £ N. We define a partition of tlie range intervai [0,R] witli nodes (?'")^^q given by r" := nk for 
n — 0, . . . , N , and let (■(i")^^Q C X/j be such that := u{r",yj) for j = 1, . . . , J — 1 and n = 0, . . . , A^, 

where u is the solution of (|1.8I) . Also, we set r"+^ := ^^+'' for n = 0, . . . , N — 1. 

Lemma 4.1. For ?j G we /lawe 

(4.5) Re(a; (X) 9^^, 'y)oji = -5 (w, Ihv)o,h- 

Proof. For w (E X^, we have 

.1-1 j-i 

.7-1 j-i .7-1 

i=i i=i j=i 
j-i 

= -(L.0 9;,t;,^)o,;.-fe^^; r^^^r^'S 

which easily yields (|4.5p . □ 
Lemma 4.2. For ?j G Xh we have 

(4.6) l|4w||o,. < \\v\\o,h- 
Proof. For w G X^, we have 



||^/iw||o,h 5 lh^\vj+i+Vj_i\' 



.7-2 



which easily yields (I4.6p . □ 
Lemma 4.3. For ?j G Xh we have 

(4.7) |i;|oo,h < V2 {\\v\\o,h + \\dhv\\o,h) ■ 

Proof. Let w G X/,, with w 7^ 0. Then, there exists io G {1, . . . , J — 1} such that |w|oo,/i = \ . If ?o is even, 
i.e. io = 2 Too for some toq G N, then we have 

bio I =l^^2mol 
rno-1 

e=o 

(4.8) 



mo — 1 \ 2 ^ mo — 1 

"2(f + l) — t'2f 
2/1 



^2 

\ £=0 / \ i=0 

<V2\\dhv\\o,h. 
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Let us now assume that io is odd, i.e. io = 2 toq — 1 for some mo € N. If J is odd, i.e. J = 2 J^, + 1 for some 
J* e N, then 

\vio \ = |w2mo-l| 



(4.9) 



<2h y: 

i—niQ 



'"2g+l~^'2f-l 

2h 



< 



2h 



i—niQ 



<V2\\dhv\\Q,h- 



Now, we assume that J is even, i.e. J = 2 J* for some J* € N. We define {w()''^*^^ C C, by wi = V2i^i for 
£ = 1, . . . , J*. Also, let e {1, . . . , J*} such that \wi^ \ = mini<f<j^ \wi\. Then, we have 



(4.10) 



<2/i ^ Iwfl + 2 ft. ^ 



2ft 



<\/2||w||o,ft + V2||9,,?;||o,ft. 

Thus, g21) follows from glU), (gH) and 

4.2. Properties of the discrete elliptic operator A?i and its inverse. 
Lemma 4.4. We assume that 



(4.11) 



Cr,EB luf 



^ (Re(g) + |gpRe(7)) 



> 0, 



_(Cdpf)2 Q2|q|2 

where C^pf G (Oj optimal constant in (|4.ip . or, t/iat fftere exists S^, G {1,-1} smc/i t/iat 



C„„„ := inf 



(Im(g)-|g|2lm(7)) 



> 0. 



max \\Thir)v\\i^h < C \\v\\o.h V« G 

rG[0,K] 



(4.12) 

Then, we have 
(4.13) 

w/iere Tft (r) := A^^(r) /or r G [0, i?]. 

Proof. Let r G [0, R] and w G AT^. Then, we have 

.7-1 

(4.14) Re(A,(r)«, v)o_, = \v\l^ - ^ h ^ [Re(g) + \q\^ Re{j{r, y,))] \v, 

When (|4TT|) holds, then (j4?T4| and (gH]) yield 



- 1 

Re(A,(r)z;,z;)o.. >hY, [(^-^ - ^ (Re(g) + |gp Re(7(r, y,))) 



(4.15) 

>Co^i.|Mlo,/,. V«GX,„ VrG[0,i?], 
When (|4.12p holds, then, we have 

I Im(A,i(r)w,w)o,,. I > (5* Im(A,i(r)w, w)o,h 



(4.16) 



> 



'^*^^EM9)-|g|'lm(7(r,j;,))] 



>CnBB\\v\\l^, yveX,,, VrG[0,i?]. 
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□ 



Now, using (j4.15|) or (|4.16p . and the Cauchy-Schwarz inequality we arrive at 

(4.17) \\Thir)v\\o,h< C\\v\\o,h. 

Next, we use (|4TT4|) and (jiTf)) to get 

\THir)v\lH < C ( \\Th{r)v\\l^ + \\v\\o,h \\Th{r)v\\o,h \ 
<C\\v\\l^. 

Thus, (j4l3l) follows easily from (j4lT| and (|4A8l) . 



(4.18) 



□ 



Proposition 4.1. We assume that (|4.1ip or (j4.12p hold. Then, there exist positive constants Ca and Cb 
such that 



(4.19) 

and 

(4.20) 



II [Th{r)Phcf,-n{T)Ph^]-[PhT{r)cl>-PhT{T)cl)]\\^,<Cs [h^ |0|3,oo,. + k - t| |,^|2,oo,. 

V<?!> e C3(/;C), Vr,Te[0,i?]. 



Then, we have that 
(4.21) 

where r]{r) € Xh with 



Proof. Let r, r G [0,i?], G C2(/;C), V'(r,y) (T(r)0)(2/) and E{r) G X,, defined by 

E{r) ■.^TH{r)PH4>-PhT{r)cp. 

Kh{r)E{r) = r]{r) 

with ^j(r') G iuj-i, Vj+i) for j = l,...,J— 1. Thus, along with (j3.10p . we obtain that 

|??(?')|oo,h < niax|?/'yyay(r,-)| 

<C/l2|<^|2,oo,.. 

We use and to have 

||i?(r)||i,„-||r;,(r)ry(r)||i,^ 
(4.23) <C||r7||o,h 

If we combine (j423l) and (|422|) . (jiTTQ]) easily follows. 
Now, assuming that r < r, we have 

('?(''))j - ('7(T))i =f? { [V'wTO(?',Ci(0) - V'TOaa(T,Cj('"))] + [^yyyv{'r,£.j{r)) - ■>l}yyyy{T,£,j{T))]} 



(4.22) 
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(r - r) a.aXA^, (r, r), (r)) + (^, (r) - (r)) ^^^^(r, (r, r)) 



, j = 1,..., J- 1, 



with /ij(r,T) G (?',t) and ^j(r,T) G (j/j_i, i/j+i) for j = 1,..., J — 1. Thus, using (|3.10p and p.lip . we obtain 
that 



(4.24) 



\r]{r) - 77(t)|ooJi < |r - t| max IV'yyaarl + 2/i Ta&y.\il}yyyyy\ 

<Ch^ (|r-T||0U,, + /i|0|3,oo,.). 
Using (|4.2ip we see that 

(4.25) Ah{r){E{r) - E{t)) = [ry(r) - r/(r)] + Z(r, r) ® Eir) 



J = 1,...,J-1. 



where Z{r, r) e is defined by 

(4.26) {Z{r,r)), :^ - [s^r] i±i^ + s^r) i±2^ 
Then, from ()4.26p we obtain 

(4.27) \Z{r,T)\^^h< C\r-T\. 
We use ffn^i and to get 

\\E{r) - E{T)h,h < C ( Mr) - r,(T)||o.„ + \\Z{r,T) ® £;(r)||o,, ) 

< C ( h(r) - vir)\oo,h + \Z{r,T)\^,n \\E{T)\\o,h ) • 
Then, follows easily from (|T?7)) and g^H). 

4.3. The Crank-Nicolson-type finite difference method CNFD. 



(4.28) 



□ 



4.3.1. Formulation of the method CNFD. For m = 0, • • • , N, the method CNFD constructs an approximation 
t/™ of li™ following the steps below: 

Step Al: Define U° € Xh by 

(4.29) [/" Phuo. 
Step A2: For n 1, . . . , iV, find C/" G such that 

(4.30) ^^^^^^ - S{r-i) u:®du (U^+u^Zl^ ^ i ^(^-i) T,(r"-^) (£::±^) + i A (^Ul+ipl' 

where <5(r) := for r e [0, R], and ^(r) := for r e [0, R]. 

4.3.2. Consistency o/ CNFD. 

Lemma 4.5. We assume that (|4.1ip or (|4.12p /lo/ii. ^Zso, /or n = 1, . . . , A^, we define rf^ G X/j hy 
""'f - (5(r"-5) ( ""+f' =i ^(r"-^)T„(r"-^) 

+ if ( ""+2"'' ) +??"• 

Then, there exists a constant C > depending on the data of the problem (|1.8|) which is such that 
(4.32) 

and 

max I77" - r,"-i|oo.^ <C' [fc^ ( \dtu\^,^ + \dluU,n + l^^^^]^ + \dldyuU^) 

2<n<N L V / 



(4.31) 



max |?7™|oo,/.< C [k" Bi{u) + h'' B2{u)] 

l<m<N 



(4.33) 



where 



+ h-^ |(9„M|oo,n + niax|u|3,oo,7 

[0,-R] 



+ kh'^ ( |i9^M|^,n + \drd^u\^^u + max|M|2,oo,/ + max|Mr|2,oo,/ 



Bi{u) :— \u. 



[0,-R] 



S2(u) := max |u|3 oo,/- 



Proof. Let n e {1, . . . , iV}. Using Taylor's formula we conclude that 



v7 = k^ 



^{T-,y,)~6{r^-^)y, 



"^yrriTj ,yj) + g Uyyy{r ,Zj)+ g Uyyylr ,Zj) 

i C(r"-^) [T,(r"-^)P„ (li{ir^±!£(lM) _ P^T{r^-^) ( "(-"v)+n(r-\-) j 
i e(r""') [PhT(r"-3) (!i(lM±|(irM _ u(r"-5, •)) 



(4.34) 



+ i^!^Urr{t],yj), i = l,...,J-l, 



where t", r", r" e (r" \r") and z" G 1,7/^+1). 



Now, combine gj!]), (g^]), (j4l9| . ([330]) and ([3lT|) to obtain 

^rrr|oo,fj |''^rr|oo,£l "t~ |'^yrr|oo,r2 j ~\- h \y'yyy\c 



Ch maxlukoo/ 

10,R] ' ' ' ' 



(4.35) 



C 



tt(r",-)+ti(r" ,-) 
2 



u(r"-2,.) 



<C 



k ( |Urrr|oo,n ~l" |^rr|oo,Q "I" |Wyrr|oo,f2 ) + ^ ^"^^ l^ls.oo,/ 



and thus arrive at (|4.32l) . 
Since it holds that 



^^7 "^1 ^ ~ 24 [ 



5{r 2) (5(r ^) yj g Uyrr{Tj ,yj)-\- g Uyyy{r , z^- ) + g Uyyy{r T^j) 



2 



for j = 1, . . . , J— 1 and n = 2, . . . , N , (|4.33p follows easily by the mean value theorem, (|4.19p and (|4.20p . □ 
4.3.3. Convergence of CNFD. 

Theorem 4.2. Let (U™)'^^q C be the finite difference approximation to the solution of (jl.8l) defined as 
in Section \4-.3.1\ Then, there exist constants Ci > 0, C2 > and C3 > independent of k and h such that: 
if C'^k < 3, then 



(4.36) 

and 

(4.37) max 

2<n<N 



max ||C/"-u"||o.h < Ci max h"|lo,h 

0<n<iv l<n<iv 



[2 2 ) 



< Co 



00. h 



fc-i max II 77" II o,h+ max h"|loj. 



Proof. Let e" it" - C/™ for n = 0, . . . , iV, and observe that due to (|4?29| there holds that e° = 0. Next, 
subtract (|430ll from (|43T|) to obtain 



(4.38) 



+ 77", n=l,...,iV. 



i A ( e"+e" 
9 I 2 
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Now, take the (•, •)o,h— inner product of both sides of (I4.38P with e" + e" ^, and then real parts to get 



(4.39) 



■Im 



f Im 



0,ft 



, n = l,...,N. 



+ fcRe(?7",e" + e"-i^ 
Combining (gSH), (g^]), (gH) and we obtain 

(4.40) |ie"|!o./.-||e"-i||o,h < kc. ( ||e"||o,ft + ||e"-i||o.O + ^ h"llo,;., n = l,...,N. 
Assuming that fc < ^ , from (I4.40p we conchide that 

nj.< ^\\e--'\kh + T^nH, n = l,...,7V, 

which yields 

(4.41) |le"||oj. < e^=*'-||e"'i||o,h + 3fc||r7"||o,h, n = l,...,iV. 

Applying a standard discrete Gronwall argument on (|4.41l) and using the fact that He^Ho,/! = 0, we obtain 



(4.42) 



max |le"||o,;,. <C(||e"||o,h + max ||77"||o,/i) 

0<n<N ^ l<ri<N " ' " ' / 



<C max ||?7"||o,/!., 

l<n<JV 



which establishes the estimate (|4.36p 
— e™"-'^ for TO = 1, . 



Let z^™ e™ - e"-i for to = 1, . . . , TV. Then, (j438l) yields 



(5(r"-5)_5(r"-2) 



(4.43) 



+ ie(^"-^)T„(r"-5)(£2±|::ii) 

+ i [e(r"-^)-e('-"-^)] T;,(r"-^) (^^^^±f^ 

+ ie(r"-^) Inir'^-i) 



-^/^(r"-*) 



(^) 



+ (7/"-77"-i), n = 2,...,iV. 
Now, take the (•, •)o,h~inner product of both sides of ()4.43p with i/" + i^"^^ and then real parts to get 



I ^i||2 II n— 1||2 ^ 1 II n I n— In , j n n , n— In , |i n 1| 



(4.44) 



T,,(r"- 2)(e" + e"~i) - Tft (r"-2 )(e" + e"-^) 



0,h 



W ( 



0,ft 



after using (|4TT3| . Now, (j438| along with (j4?T3| . yield 



< C 



0,h 



n = 2,...,iV. 



(4.45) juj<S)dh ( ''"'t""'' ) 
Thus, (|4:45l) and (|4:44l) yield 

lk"llo,h - lk""'||o./. <Ck[ ||i."||oj. + |k"-'||oj. + k {\\e"\\oM + \\e"-%,h + ||e""'||o,/.) 

(4.46) +A:||ry"-i||o.ft + h"-r7"-i||oj. 



+ 



Th(r"-^)(e" + e"-^) - T;,(r"-5)(e« + e""^) 
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for n = 2, . . . , TV. Now, we observe that, for v e Xh, we have 
(4.47) ^ ^ 

= V^"(^r,,(r"-i)i;, n = 2,...,7V, 

where V" € is given by 



(4.48) V- :^ .^(r"-i) - s^(r""^) ^^^'^l:""" j = l,...,J-l. 



Using (|iTf|) . and we have 

<Ck\\v\\o,h, yveXh, n = 2,...,iV. 

Combining (|4.46p and (|4.49p we obtain 



(4.50) 



lk"||o,h - lk"^'llo,/. < a fc |k"||oj. + \W"-'\\o,k + k {n,k + \\e"-%M + \K-^\\om) 



+ fch""'l|o,/, + h"-??"-'||oj. 



n = 2, . . . , A^. 



Assuming that k is enough smah (i.e. 3fc max{c*,C*} < 1) and applying a discrete Gronwah argument on 
(|4.50p . we conclude that 



max \\iy"-\\o.h < C 

l<n<N 



max |l?7"-?7" ^\\o,h + k max |l?7"|lo./i 

2<n<Af l<n<Af-l 



+ k max ||e"||o,/i + ||j^^||o,/i 

0<n<N . 



which, along with (|4.36p and (|4.4ip . yields 



(4.51) 



max ||i/"||o,h<C max h" - 77""^ llo./i + fc max h"|lo,/i + |le^|lo,/i 

l<n<iv L2<n<iv l<n<iv 



<c 



max II?/" -r/" -^||o,/i + fc max ||?y"||o,/i 

2<n<N l<n<N 



Now, from (|438)) follows that 
(4.52) ||c^ ® dh{e'' + e"-i)||o,„ < C [A:"! ||i."||o,,, + ||e"||o,,, + ||e"-^||o,/. + ||?/"||o,/. 
Then, (|432)) and (|43Tt yield that 



n = l,...,Af. 



(4.53) max ||w 9,,(e" + e"-^)||o,,, < C A;-^ max - t?""' ||o,h + max ||77"||o,,, 

l<n<N L 2<71<N 1<71<N 



Finally, (|4.37p follows from (14. 53^ and (|4.7p . in view of the identity dh{i^<E)v) — dhV + IhV for v G AT/j. □ 

In view of the results of Lemma 14. 5[ (|4.36p implies that 

max ||[/"-u"||o, = C'(/c2 + /i2). 

0<ra<JV 

Also, the estimates (|i37|) . and yield 

max |w(g)e"-3|^^ = Oik^ + + k-^) 

0<ri<JV 

which may be viewed as a discrete weighted maximum norm estimate on e"~ ^ , which is of optimal order 
when h = 0{k); however does not yield an optimal order maximum norm estimate, due to vanishing of the 
coefficient of Uy in the p.d.e. in ()1.8p at y = 0. 
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5. Numerical implementation 
5.1. The numerical scheme. Using the definitions S{r) :— CI''; 2/) ■— ^ '^''^ ^ ^(y) :— 



we note that the p.d.e. in ()1.8p may written as 

-C('', y) - dyG{r, y) = if (r) u, 

where we have put G(r, y) := Ur{r, y) — i — u{r, y) — 5{r)y Uy{r, y). Motivated by this we rewrite the CNFD 
scheme (I4.29p - (l4.30p in the fohowing equivalent form using the notation introduced in section 2] and putting 

^ :— — — 2~— We seek [/", < j < J, < n < TV, approximating u" and given by the equations: 
For n = 0: 

(5.1) U^ = Uoiyj) for 1<J<J-1, C/0 = C/," = 0. 
For n = 1, . . . ,iV: 

(5.2) -C(r"-iy,)G,"-A„G^"=.ie(r"-^)f/;"^ for I < j < J ~ 1, 
where 



/ G£-2G£ 

—fP- — 

-2G'; ,+G'; , 



i = 1, 

J = 2, • • • , J - 2, 



with 



Gn,^Hlz^-ilu;-^-5{r--^)y,^^^±^^P^ for j = l,...,J-l. 
The scheme requires solving a pentadiagonal linear system of algebraic equations at each time step. 



5.2. Numerical experiments. We implemented the finite difference scheme CNFD in the form (|5.ip - (|5.2p 
in a double precision FORTRAN 77 code using it in various numerical examples to test its accuracy and 
stability. We used the function u{r, y) = exp(2r)(y— 1) sin(27ry) as exact solution of (II. 8p (with an appropriate 
nonhomogeneous term in the right-hand side of the p.d.e.), putting 7(r, y) ==l + j/, q: = 2, p — q+i, and 
selecting q = (0.252252311, —1.35135138 e — 2), [T^. We experimented with several bottom profiles. For 
example, in the case of the downsloping bottom given by s(r) = exp(r), < r < 1, we obtained at r = 1 
the errors in the discrete and norms shown in Table [T] together with the associated experimental 
convergence rates. (We took h = k = ^ for the values of J shown). The convergence rates of the table are 
practically equal to 2 and consistent with the predictions of the theory. We also tried bottom profiles s(r) 
given for < r < 1 by r + 2, — r + 2, — exp(— r), cos(27rr) + 2, i.e. upsloping and oscillating profiles as 
well, and found experimentally that the scheme was stable ior k — h and that the L^- and convergence 
rates were practically equal to 2 again. Hence it seems that although CNFD is designed to approximate the 
ibvp ()1.8p in the downsloping bottom case, it is resilient enough in upsloping and non-monotonic bottom 
problems as well. 

Although the ibvp (II. 8p is L^-conservative, in the sense that its solution preserves the quantity 

\/s{?j\\u{r,-)\\mi) 



J 


L'^-crror 


L'^-rate 


L°° -error 


L°°-rate 


40 


0.2510e-l 




0.2493 e-1 




80 


0.6424 e-2 


1.966 


0.6365 e-2 


1.969 


160 


0.1627e-2 


1.981 


0.1609 e-2 


1.983 


320 


0.4097 e-3 


1.990 


0.4048 e-3 


1.991 


640 


0.1028e-3 


1.995 


0.1015e-3 


1.995 


1280 


0.2574 e-4 


1.997 


0.2542 e-4 


1.998 



Table 1. Discrete L^- and L°°-errors and rates at r = 1, s(r) — exp(r), k = h = j. 
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for real 7 and q, CNFD does not share a corresponding discrete property. For example, when we integrated 
numerically (jl.8|) with a — 10, q — j, p — q + ^, 7(r, y) = 1 + y, uo{y) = y'^(y — 1), < y < 1, we found 
that in the cases s{r) = exp(r) and s(r) = r + 2, 0<7'<1, the quantity Vs(r^||[/"||o,h was preserved to 4 
significant digits. 

We also performed a simulation of a realistic underwater acoustics problem using the CNFD scheme. We 
integrated the ibvp (|1.8p using again q — (0.252252311, —1.35135138 e — 2), p — q + ^, and considered a 
straight downsloping bottom given by s{r) = 200 (l + 4500) (distances in meters) and making an angle of 
2.86° with respect to the horizontal surface. We used as initial condition at r = the normal mode starter 
given by formula (45) of |12| with M = 6, simulating the initial field produced by a time-harmonic point 
source of frequency f = 25 Hz located at a depth of 2: = 100 m. (We shall frequently refer in the sequel 
to quantities expressed in the r, z variables. Of course the scheme is implemented in the r, y variables 
and the results are transformed from or into the r, z domain as required.) We assume that the medium is 
homogeneous and lossless with a sound speed c = Cq = 1500 m/sec, so that P = 0. We integrated the problem 
up to i? = 3300 m using k — 0.83475 m, and 4000 mesh intervals of equal length in y. As is customary in 
underwater acoustics we present the numerical results in terms of a one-dimensional transmission loss (TL) 

plot in the r, z variables. (The TL function was computed by the formula TL= —20 logj^g ^ ^ where 

Zrec is a receiver depth.) The graph of Figure [2] shows as an example, the TL curves at Zroc = 30 m obtained 
by integrating CNFD (solid line) and by the finite difference scheme (dotted line, 'DSZ scheme') proposed in 
[12) as a discretization of a problem of the form (|1.8|) but with the bottom boundary condition (cf. Remark 
[33D 

Uyy{r,l) =if s(r)s(r)wy(r, 1), < r < i?, 
replacing Uy{r, 1) — 0. The results of the two schemes are in very good agreement. (In T^, the results of 
the DSZ scheme were compared to those of a standard 'staircase' wide-angle PE code and were found to be 
in very good agreement.) 
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Figure 2. TL at z^-cc = 30 m, downsloping bottom. CNFD (solid line), DSZ (dotted line). 
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